% 加载数据
data = load('my_signal_with_noise5.mat');
fs = data.fs;  % 获取采样频率
signal_noisy = data.signal_noisy;  % 获取带噪声信号  

% 选择 PCA 或 ICA 方法提取信号
% 假设我们已经使用 PCA 或 ICA 恢复了信号，这里以 PCA 为例
window_size = 100;  % 将信号分为多个窗口
X = buffer(signal_noisy, window_size, window_size - 1, 'nodelay');  % 将信号分割为窗口
[coeff, score, ~] = pca(X');  % 使用 PCA 提取主成分
reconstructed_signal = score(:, 1) * coeff(:, 1)';  % 重建信号
%reconstructed_signal = reconstructed_signal(:);  % 还原为列向量

% 计算恢复后信号的频率和相位信息
N = length(reconstructed_signal);
f = (0:N-1) * (fs / N);  % 频率轴（从0到采样频率）

% 快速傅里叶变换（FFT）
signal_fft = fft(reconstructed_signal);
magnitude = abs(signal_fft);  % 幅度
phase = angle(signal_fft);    % 相位

% 只取前半部分的频谱（因为FFT是对称的）
f = f(1:N/2);
magnitude = magnitude(1:N/2);
phase = phase(1:N/2);

% 找出最大幅度对应的频率（主要频率成分）
[peak_magnitude, peak_index] = max(magnitude);
peak_frequency = f(peak_index);  % 主频率
peak_phase = phase(peak_index);  % 对应的相位

% 绘制恢复后的信号及其频率和相位
t = (0:N-1) / fs;

signal = data.signal;
% 绘制时域信号
subplot(2, 1, 1);
plot(t, reconstructed_signal,'b');
title('Denoised Signal');
xlabel('Time [s]');
ylabel('Amplitude');


% % 绘制频率谱（幅度谱）
% subplot(4, 1, 2);
% plot(f, magnitude);
% title('Frequency Spectrum');
% xlabel('Frequency [Hz]');
% ylabel('Magnitude');
% 
% % 绘制相位谱
% subplot(4, 1, 3);
% plot(f, phase);
% title('Phase Spectrum');
% xlabel('Frequency [Hz]');
% ylabel('Phase [radians]');


N2=length(signal);
t2 = (0:N2-1) / fs;
subplot(2, 1, 2);
plot(t2, signal,'b');
title(' Signal');
xlabel('Time [s]');
ylabel('ll');



% 输出主频率和相位信息
fprintf('Peak Frequency: %.2f Hz\n', peak_frequency);
fprintf('Peak Phase: %.2f radians\n', peak_phase);
